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Abstract 

Construction of an accurate theory of orbits about a precessing and nutating oblate 
planet, in terms of osculating elements defined in a frame associated with the equator of 
date, was started in Efroimsky & Goldreich (2004) and Efroimsky (2004, 2005, 2006a,b). 
Here we continue this line of research by combining that analytical machinery with numerical 
tools. Our model includes three factors: the J2 of the planet, its nonuniform equinoctial 
precession described by the Colombo formalism, and the gravitational pull of the Sun. This 
semianalytical and seminumerical theory, based on the Lagrange planetary equations for the 
Keplerian elements, is then applied to Deimos on very long time scales (up to 1 billion of 
years). In parallel with the said semianalytical theory for the Keplerian elements defined 
in the co-precessing equatorial frame, we have also carried out a completely independent, 
purely numerical, integration in a quasi-inertial Cartesian frame. The results agree to within 
fractions of a percent, thus demonstrating the applicability of our semianalytical model over 
long timescales. 



* We use the term "precession" in its general meaning, which includes any change of the instantaneous spin 
axis. So generally defined precession embraces the entire spectrum of spin-axis variations - from the polar wander 
and nutations through the Chandler wobble through the equinoctial precession. 
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Another goal of this work was to make an independent check of whether the equinoctial- 
precession variations predicted for a rigid Mars by the Colombo model could have been 
sufficient to repel its moons away from the equator. An answer to this question, in combi- 
nation with our knowledge of the current position of Phobos and Deimos, will help us to 
understand whether the Martian obliquity could have undergone the large changes ensuing 
from the said model (Ward 1973; Touma & Wisdom 1993, 1994; Laskar & Robutel 1993), or 
whether the changes ought to have been less intensive (Bills 2006, Paige et al. 2007). It has 
turned out that, for low initial inclinations, the orbit inclination reckoned from the precess- 
ing equator of date is subject only to small variations. This is an extension, to non-uniform 
equinoctial precession given by the Colombo model, of an old result obtained by Goldreich 
(1965) for the case of uniform precession and a low initial inclination. However, near-polar 
initial inclinations may exhibit considerable variations for up to ±10 deg in magnitude. 
This result is accentuated when the obliquity is large. Nevertheless, the analysis confirms 
that an oblate planet can, indeed, afford large variations of the equinoctial precession over 
hundreds of millions of years, without repelling its near-equatorial satellites away from the 
equator of date: the satellite inclination oscillates but does not show a secular increase. Nor 
does it show secular decrease, a fact that is relevant to the discussion of the possibility of 
high-inclination capture of Phobos and Deimos. 



1 Introduction 

1.1 Statement of purpose 

The goal of this paper is to explore, by two very different methods, inclination variations of 
a solar-gravity-perturbed satellite orbiting an oblate planet subject to nonuniform equinoctial 
precession. This nonuniformity of precession is caused by the presence of the other planets. Their 
gravitational pull entails precession of the circumsolar orbit of our planet; this entails variations 
of the solar torque acting on it; these torque variations make the planet's equinoctial precession 
nonuniform; and this nonuniformity, in its turn, influences the behaviour of the planet's satellites. 
This influence is feeble, and we trace with a high accuracy whether it results, over aeons, in purely 
periodic changes in inclination or can accumulate to secular changes. 

This work is but a small part of a larger project whose eventual goal is to build up a compre- 
hensive tool for computation of long-term orbital evolution of satellites. Building this tool, block 
by block, we are beginning with only three components - the planet's oblateness, the direct pull 
of the Sun on the satellite and the planet's precession. These phenomena bare a marked effect 
on the evolution of the orbit. In our subsequent publications, we shall incorporate more effects 
into our model - the triaxiality, and the bodily tides. 

1.2 Motivation 

One motivation for this work stems from our intention to carry out an independent check of 
whether the equinoctial-precession changes predicted for a rigid Mars by the Colombo model 
could have been sufficient to repel its moons away from the equator. An answer to this question, 
in combination with our knowledge of the current position of Phobos and Deimos, will help us 
to understand whether the Martian obliquity variations could indeed have undergone the large 
variations resulting from the Colombo model, or whether the actual variations ought to have 
smaller magnitude. Such a check is desirable because the current, Colombo-model-based theory 
of equinoctial precession (Ward 1973; Touma & Wisdom 1993, 1994; Laskar & Robutel 1993), 
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incorporates several approximations. First, the Colombo equation is derived under the assertion 
that the planet is rigid and is always in its principal spin state, the angular-momentum vector 
staying parallel to the angular-velocity one. Second, this description, being only a model, ignores 
the possibility of planetary catastrophes that might have altered the planet's spin mode. Third, 
this description ignores that sometimes even weak dissipation (caused, for example, by tides) 
may be sufficient to quell chaos and regularise the motion, which may be the case of Mars (Bills 
2006). Fourth, it still remains a matter of controversy as to whether the observed pattern of small 
craters on Mars confirms (Hartmann 2007) or disproves (Paige et al. 2007) the strong climatic 
variations predicted in Ward (1974, 1979). All this motivates us to come up with a test based 
on the necessity to reconcile the variations of spin with the present near-equatorial positions of 
Phobos and Deimos. The fact that both moons found themselves on near-equatorial orbits, in all 
likelihood, billions of years ago,^ and that both are currently located within less than 2 degrees 
from the equator, is surely more than a mere coincidence. An elegant but sketchy calculation by 
Goldreich (1965) demonstrated that the orbits of initially near-equatorial satellites remain close 
to the equator of date for as long as some simplifying assumptions remain valid. As explained in 
Efroimsky (2004, 2005), these assumptions are valid over time scales not exceeding 100 million 
years, while at longer times a more careful analysis is required. Its goal will be to explore the 
limits for the possible secular drift of the satellite orbits away from the evolving equator of date. 
Through comparison of these limits with the present location of the Martian satellites, we shall be 
able to impose restrictions upon the long-term spin variations of Mars. If, however, it turns out 
that near-equatorial satellites can, in the face of large equinoctial-precession variations, remain 
for billions of years close to the moving equator of date, then we shall admit that Mars' equator 
could indeed have precessed through billions of years in the manner predicted by the Colombo 
approximation. 

The second motivation for our study comes from the ongoing discussion of whether the Martian 
satellites might have been captured at high inclinations, their orbits having gradually approached 
the equator afterwards. While a comprehensive check of this hypothesis will need a more detailed 
model - one that will include Mars' triaxiality, the tidal forces, (Lainey, Gurfil, & Efroimsky 
2008), and perhaps other perturbations - the first, rough sketch of this test can be carried out 
with only J2 , the Sun, and the equinoctial precession taken into account. We perform such a 
rough check for a hypothetical satellite that has all the parameters of Deimos, except that its 
initial inclination is 89° . 

^ Phobos and Deimos give every appearance of being captured asteroids of the carbonaceous chondritic type, 
with cratered surfaces older than ~ 10^ years (Veverka 1977; Pang et al. 1978; Pollack et al. 1979; Tolson et al. 
1978). If they were captured by gas drag (Burns 1972, 1978), this must have occurred early in the history of the 
solar system while the gas disk was substantial enough. Kilgorc, Burns, and Pollack (1978) have demonstrated 
numerically that a gas envelope extending to about ten Martian radii, with a density of 5 x 10~^g/cm'^ at the 
Martian surface, could have been capable of capturing satellites of radii about 10km. At that stage of planetary 
formation, the spin of the forming planet would be perpendicular to the planet's orbit about the Sun ~ i.e., the 
obliquity would be small and the gas disk would be nearly coplanar with the planetary orbit. Energetically, a 
capture would likely be equatorial. This is most easily seen in the context of the restricted three-body problem. 
The surfaces of zero velocity constrain any reasonable capture to occur from directions near the inner and outer 
coUincar Lagrange points (Szebehely 1967, Murison 1988), which lie in the equatorial plane. Also, a somewhat 
inclined capture would quickly be equatorialised by the gas disk. If the capture inclination is too high, the orbital 
energy is then too high to allow a long enough temporary capture, and the object would hence not encounter 
enough drag over a long enough time to effect a permanent capture (Murison 1988). Thus, Phobos and Deimos 
were likely (in as much as we can even use that term) to have been captured into near-equatorial orbits. 
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1.3 Mathematical tools 

The first steps toward the analytical theory of orbits about a precessing and nutating Earth 
were undertaken almost half a century ago by Brouwer (1959), Proskurin & Batrakov (1960), 
and Kozai (1960). The problem was considered, in application to the Martian satellites, by 
Goldreich (1965) and, in regard to a circumlunar orbiter, by Brumberg, Evdokimova, & Kochina 
(1971). The latter two publications addressed the dynamics as seen in a non-inertial frame co- 
precessing with the planet's equator of date. The analysis was carried out in terms of the so-called 
"contact" Kepler elements, i.e., in terms of the Kepler elements satisfying a condition different 
from that of osculation. Modeling of perturbed trajectories by sequences of instantaneous ellipses 
(or hyperbolae) parameterised with such elements is sometimes very convenient mathematically 
(Efroimsky 2006c). However, the physical interpretation of such solutions is problematic, because 
instantaneous conies defined by nonosculating elements are nontangent to the trajectory. Though 
over restricted time scales the secular parts of the contact elements may well approximate the 
secular parts of their osculating counterparts (Efroimsky 2004, 2005), the cleavage between them 
may grow at longer time scales. This is the reason why a practically applicable treatment of the 
problem must be performed in the language of osculating variables. 

1.4 The plan 

The analytical theory of orbits about a precessing oblate primary, in terms of the Kepler elements 
defined in a co-precessing (i.e., related to the equator of date) frame, was formulated in Efroimsky 
(2006a,b) where the planetary equations were approximated by neglecting some high-order terms 
and averaging the others. This way, from the exact equations for osculating elements, approximate 
equations for their secular parts were obtained. We shall borrow those averaged Lagrange-type 
planetary equations, shall incorporate into them the pull of the Sun, and shall numerically explore 
their solutions. This will give us a method that will be semianalytical and seminumerical. We 
shall then apply it to a particular setting - evolution of a Martian satellite and its reaction to the 
long-term variations of the spin state of Mars. Our goal will be to explore whether the spin-axis 
variations predicted for a rigid Mars permit its satellites to remain close to the equator of date 
for hundreds of millions through billions of years. In case the answer to this question turns out 
to be negative, it will compel us to seek nonrigidity-caused restrictions upon the spin variations. 
Otherwise, the calculations of the rigid-Mars inclination variations will remain in force (and so 
will the subsequent calculations of Mars obliquity variations); this way, the theory of Ward (1973, 
1974, 1979, 1982), Touma & Wisdom (1993, 1994), and Laskar & Robutel (1993) will get a 
model-independent confirmation. 

2 Semianalytical treatment of the problem 

To understand the evolution of a satellite orbit about a precessing planet, it is natural to model it 
with elements defined in a coordinate system associated with the equator of date, i.e., in a frame 
co-precessing (but not co-rotating) with the planet. A transition from an inertial frame to the 
co-precessing one is a perturbation that depends not only upon the instantaneous position but 
also upon the instantaneous velocity of the satellite. It has been demonstrated by Efroimsky & 
Goldreich (2004) that such perturbations enter the planetary equations in a nontrivial way: not 
only they alter the disturbing function (which is the negative Hamiltonian perturbation) but also 
they endow the equations with several extra terms that are not parts of the disturbing function. 
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Some of these nontrivial terms are linear in the planet's precession rate /x , some are quadratic 
in it; the rest are linear in its time derivative /* . The inertial-forces-caused addition to the 
disturbing function (i.e., to the negative Hamiltonian perturbation) consists of a term linear and 
a term quadratic in fi . (See formulae (53 - 54) in Efroimsky (2004, 2005) or formulae (1) and 
(6) in Efroimsky (2006a,b).) The essence of approximation elaborated in Ibid, was to neglect 
the quadratic terms and to substitute the terms linear in /x and /x with their secular parts 
calculated with precision up to , inclusively. 

2.1 Equations for the secular parts of osculating elements defined in 
a co-precessing reference frame. 

We shall begin with five Lagrange-type planetary equations for the secular parts of the orbital 
elements defined in a frame co-precessing with the equator of date. These equations, derived in 
Efroimsky (2004, 2005), have the following form: 

? = - 2^a , (1) 

dt n ^ ^ ' ^ ^ 

^ = ^^e (1 - e^Y^' , (2) 
dt 2 n ^ ^ ^ ' 

duj 3 nJ2 fPe\^f5 2 ■ 1\ , cosz /■ f f^^^\\ ro\ 

- ' ' - cos I - - ] - n± + Hn coil J- ^ . — : (/X -f X — ) , (3) 



dt 2(l-e2)2Vay V2 2j ' ^ na^ {1 - e^^^ sim \ di 



— = — Ml cosVL — ji2 sinf2 
dt 



+ 



cos I 



na^(l — e^)-*^/^ sin i 



[1 - e2)V2 sin 



(4) 



^ _ _ 3 J f PeV COSZ /i^ 1 /■[_{■ \ 

dt ~ 2""^ \ a J [1 - e^f sin? na2 (1 - e2)i/2 sin^ V di J ^ ' 

The number of equations is five, because one element. Mo , was excluded by averaging of the 
Hamiltonian perturbation and of the inertial terms emerging in the right-hand sides. In the 
equations, n = ^/ G {mpHmary + msecondary) /cfi , while f(t; a, e, u, Q, Mo) is the 
implicit function that expresses the unperturbed two-body dependence of the position upon the 
time and Keplerian elements. Vector /x denotes the total precession rate of the planetary equator 
(including all spin variations - from the polar wander and nutations through the Chandler wobble 
through the equinoctial precession through the longest-scale spin variations caused by the other 
planets' pull), while fii , fj,2 , Ats stand for the components of /x in a co-precessing coordinate 
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system x , y , 2; , the axes x and y belonging to the equator-of-date plane, and the longitude of 
the node, Q , being reckoned from x : 

/X = /ii X + /i2 y + /is z , where /ii = Jp , /i2 = hp sin Ip , /is = hp cos Ip , (6) 

Jp , hp being the inclination and the longitude of the node of the equator of date relative to 
that of epoch, and a dot standing for a time derivative. The quantity fi± is a component of 
H directed along the instantaneous orbital momentum of the satellite, i.e., perpendicular to the 
instantaneous osculating Keplerian ellipse. This component is expressed with 

fj,± = /Li ■ w = /ii sin2 sinl] — /i2 sin i cosfl + /is cos i , (7) 

the unit vector 

w = x sin i sin Q — y sin i cos + z cos i (8) 

standing for the unit normal to the instantaneous osculating ellipse. Be mindful that fi± is 
defined not as d{iJ, ■w)/dt but as 

/i_L = fi w = fli sin i sin Vl — ^2 sin i cos + /is cos i . (9) 

The quantity /i„ is a component pointing within the satellite's orbital plane, in a direction 
orthogonal to the line of nodes of the satellite orbit relative to the equator of date: 

l^n = ~ fJ'i sin Q cos ^ + /i2 cos Q cos i + /is sin i 

(10) 

= — Ip sin Q cos i + hp sin Ip cos Q cos i + hp cos Ip sin i . 
Its time derivative taken in the frame of reference co-precessing with the equator of date is: 
An = — /ii sin VL cos i + fi2 cos Q cos i + /is sin i 

= — ip sin fl cos i + (^ hp sin Ip + hp Ip cos Ip j cos fl cos i + (^ 'hp cos Ip — hp Ip sin Ip j sin i . 

As shown in Efroimsky (2006a,b), the /t-dependent terms, emerging in equations (1 - 5), are 
expressed with 



— I /ii [ — (2 + 3e^) cosfl + 5e^ (cosfi cos2ct; — sinfi sin2ci; cos «) ] + 



/i2 [ — (2 + 3e^) sinf2 + 5e^ (sinf2 cos2co' + cosfi sin2co' cos ^) ] + 



/is [ 5 sin2ci; sin z ] } (12) 



;ii) 
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( /t ■ I — f X — I ) = — — (2 + 3 e^) ( /ii sin z sin f2 — fi2 sin i cos + fis cos i ) 



a2 



^ { /ii sin z [ — (2 + 3 e^) sin f2 cos i + 5 ( cos ^2 sin 2uj + sin cos 2a; cos i ) ] 



+ /i2 sin i [ (2 + 3e^) cosfi cos« + 5e^ (sinfi sin2u; — cosf2 cos 2a; cosz) ] 



- /is [ (2 + 3 6^) (2 - sin^ z) + Se^ sin^ % cos 2a; ] } ( 

To integrate equations (1 - 5), with expressions (12 - 14) inserted therein, we shall need to know, 
at each step of integration, the components of fx and jj, . 

2.2 Calculation of the components of fi and fi . 

At each step of our integration, the components of ^ and ft will be calculated in the Colombo 
approximation. Physically, the essence of this approximation is two- fold: first, the solar torque 
acting on the planet is replaced by its average over the year; and, second, the precessing planet 
is assumed to be always in its principal spin state. While a detailed development (based on the 
work by Colombo (1966)) may be found in the Appendix to Efroimsky (2006a,b), here we shall 
provide a concise list of resulting formulae to be used. 

The components of fj, are connected, through the medium of (6), with the inclination and the 
longitude of the node of the moving planetary equator, Ip and hp , relative to some equator of 
epoch. These quantities and their time derivatives are connected with the unit vector k aimed 
in the direction of the major-inertia axis of the planet: 

k = ( sin/p sin hp , — sin Jp cos hp , cos/p ) . (15) 
This unit vector and its time derivative 



dk 

It 



— = ( Jp cos Ip sin hp + hp sin Ip cos hp , — Ip cos Ip cos hp + hp sin Ip sin hp , — Ip sin Ip 



depend, through the Colombo equation 

dk 



a (n-k) (kxn) , (17) 



upon the unit normal to the planetary orbit. 



T 

n = (sin lorb sin florb , — sin lorb cos florb , COS lorb ) 
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Qorb and lorb being the node and inclination of the orbit relative to some fiducial fixed plane, 
and a being a parameter proportional to the oblateness factor J2 . In our computations, we 
employed the present-day^ value of a - see Table 5 below. We chose this value because it was 
the one used by Ward (1973, 1974), and we wanted to make sure that our plot for obliquity 
evolution. Fig. 3, coincided with that of Ward (1974). 

To find the components of /x , one must know the time evolution of hp and Ip , which can 
be determined by solving a system of three differential equations (17), with VLorb and lorb being 
some known functions of time. These functions may be computed via the auxiliary variables 

q = sin lorb sinQorb , P = sin lorb cosQorb , (19) 

whose evolution will be given by the formulae: 

00 

q = Y,^J sin {s'jt + 6j) , (20) 
i=i 

00 

p = 5^ iV, cos [s'^t + 6,) . (21) 
i=i 

The following choice of the values of the amplitudes, frequencies, and phases will make equations 
(19 - 21) render Qorb and lorb relative to the solar system's invariable plane (Ward, 1974): ^ 



J 


N, 


s'j [arcsec/yr] 


S'i [deg] 




1 


0.0018011 


-5.201537 


272.06 


2 


0.0018012 


-6.570802 


210.06 


3 


-0.0358910 


-18.743586 


147.39 


4 


0.0502516 


-17.633305 


188.92 


5 


0.0096481 


-25.733549 


19.58 


6 


-0.0012561 


-2.902663 


207.48 


7 


-0.0012286 


-0.677522 


95.01 



Table 1: Numerical values used in Ward's model of the inclination and node of the Martian orbit 

The development by Ward (1974) is limited in precision. A more accurate treatment was 
offered by Laskar (1988). At the future stages of our project, when developing a detailed physical 
model of the satellite motion, we shall employ Laskar's results. In the current paper, we are 
checking if the orbital averaging of the precession-caused terms is permissible at large time scales. 
Hence, for the purpose of this check, as well for the qualitative estimate of the long-term behaviour 
of the satellites, we need a realistic, not necessarily highly accurate, scenario of the precession 
variations. 

^ An accurate treatment shows that due to the precession of the Martian orbit a exhibits quasi-periodic 
variations of about 3 % over long time scales. Neglecting this detail in the current work, we shall take it into 
account at the further stage of our project (Lainey, Gurfil & Efroimsky 2008). 

Ward (1974) has calculated the values of the coefficients TV, s', S based on the work of Brouwer & van 
Woerkom (1950) who had calculated the values of p and q relative to the ecliptic plane of 1950. Ward (1974) 
transformed Brouwer & van Woerkom's values to the solar system's invariable plane. Be mindful that, no matter 
what the reference plane, both Ward (1974) and Brouwer & van Woerkom (1950) used the same epoch, J1950. 
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2.3 The Goldreich approximation 



The above semianalytical treatment not only yields plots of the time dependence of the mean 
elements but also serves as a launching pad for analytical approximations. For example, an 
assumption of a and e being constant, and a neglect of the /t-dependent terms in (4-5), as 
well as of the term sin i in (5), gives birth to the Goldreich (1965) approximation: 



da 

'di 

de 

dt 
di 

di 

dVL 3 /pe\2 cosz 

^"^2 - 7^ , (25) 



, (22) 
, (23) 
— /ii cos Q — fi2 sin Q , (24) 



dt 2 ^ \ aJ (1 - e2)2 

du; 3 fpe\'^ 5 cos^z - 1 u„ cos z 

7 n J2 — — + — — /i± , (26) 



dt 4 V a / (1 ~ 6^)2 sin i 

the equinoctial precession being assumed uniform: 

ip = (27) 

hp = — a cos Ip (28) 

fti = (29) 

/i2 = hp sin Ip (30) 

/i3 = /ip cos/p . (31) 

For the details of Goldreich's approximation see also subsection 3.3 in Efroimsky (2005). 



2.4 The Gravitational Pull of the Sun 

Reaction of a satellite on the planetary-equator precession is, in a way, an indirect reaction of the 
satellite on the presence of the Sun and the other planets. Indeed, the pull of the other planets 
makes the orbit of our planet precess, which entails variations in the Sun-produced gravitational 
torque acting on the planet. These variations of the torque, in their turn, lead to the variable 
equinoctial precession of the equator, precession "felt" by the satellite. It would be unphysical 
to consider this, indirect effect of the Sun and the planets upon the satellite, without taking into 
account their direct gravitational pull. In this subsection, we shall take into account the pull of 
the Sun, which greatly dominates that of the planets other than the primary. 

In what follows, m and m' will be the masses of the satellite and the Sun, correspondingly, 
r and r' will stand for the planetocentric positions of the satellite and the Sun, S will signify 
the angle between these vectors. Then the Sun-caused perturbing potential Rs , acting on the 
satellite, will assume the form of 

Rs = Gm'(^ ' -'-^] (32) 

This can be expanded in a usual manner over the Legendre polynomials of the first kind. Since 
r' ^ r , we shall take only the first term in the series:^ 

^ This is justified since the next term in the Legendre series is about 2 orders of magnitude smaller than the 
potential variations generated by the precession. 
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Rs 



Gm' 



^Pl(cOs5)+ (^)'P2(C0S5) 



r cos S 



(33) 



As the term Gm' /r' is not dependent of the satelhte's elements, one has only to consider the 
restricted potential 



Gm' 



(^)'p,(cos^) 



n'2a2 fa'^^ 



— (3 cos' S 



(34) 



n' and a' being the mean motion and the semi-major axis of the Sun. 

To obtain the Lagrange-type equations for the third-body perturbation, we must first derive 
an expression for the angle S. To that end, define the directional cosines ^ = P ■ f ' and 
6 = Q ■ f ' , where r ' is a unit vector pointing from the planet toward the Sun, while P and 
Q are unit vectors of a perifocal coordinate system associated with the osculating orbital plane 
of the satellite, with P pointing to the instantaneous periapse and Q being orthogonal to P. 
Assuming that the planet's orbit about the Sun is circular, we arrive at 

^ = cosu cos{Q — M') — cost sinu sin{Q — M') , (35) 

6 = — sinu cos(f2 — M') — cos i cosu; sin(r2 — M') , (36) 

where M' is the mean anomaly of the Sun in the planetocentric frame. With aid of these relations, 
the angle S may be written down as 

cos S = ^ cos u + 9 smu , (37) 

u being the true anomaly of the satellite in the planetocentric coordinate system. Substituting 
(35) and (36) into (37), and averaging over the satellite's mean anomaly, we arrive at the Lagrange- 
type equations: 



da 

de 
d^ 

di 
d^ 

dio 
dr 







2 sin i 



(38a) 



sin' i sin 2u; + (2 — sin' i) sin 2 uj cos 2Vt + 2 cos i cos 2uj sin 2^2 (38b) 



2 



5 e' cos ^ sin2u; (1 — cos 217) — [2 + e'(3 + 5 cos 2a;)] sin 2^2 > (38c 



4 + e — 5 sin z + 5 (sin i 



cos2cij + 5 (e' — 2) cos z sm2uj sin 2^2 



+ [5(2 - e' - sin'z) cos 2u; - 2 - 3 e' + 5 sin' z] cos21]| 



(38d) 



dVt 

dr 



K — 



L= |[2 + e'(3-5 cos2cc')] (1- cos2i7) cosii-5e' sin2u; sin21]| (38e) 



Vl-e 

where we used the following set of notations: 

r = P n {t — to) 



' ^3 



P 



3m a 
16 Ma' 



(39) 
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16n / M\ 
Sn' \ m! J 



(40) 



= - A' 



(41) 



A' being the mean longitude of the Sun in the planet's frame, and M being the mass of Mars. 

An additional averaging can be performed over the motion of the planet about the Sun. 
Mathematically, this is the same as averaging over the motion of the Sun about the planet - in 
both cases averaging over A' is implied. This averaging (Innanen et al. 1997) will simplify (38) 
into 



— = 10 e Vl — /? n sin^ i sin luj 



da 

de 

di 

di 

di 
duj 
H 
dQ 







10 e2 



sm I 



2 



(3 n cos i sin 2uj 



VT 



VT 



/? n [4 + — 5 sin^ i + 5 (sin^ i — e^) 



P n [ 2 + (3 — 5 cos 2u) ] cos i 



cos 



2u] 



(42a) 
(42b) 
(42c) 
(42d) 
(42e) 



It is important to note that the calculations leading to Eqs. (38), and to their double-averaged 
version, Eqs. (42), were performed in the Martian-orbital, and not Martian-equatorial frame, 
without taking either the Martian obliquity or precession into consideration. Had we taken into 
account the precession, we would get, on the right-hand side of (42) resonances between the 
motion of the Sun relative to the planet and the equinoctial precession. Since the time scale of 
the former exceeds, by orders of magnitude, the time scale of the latter, we may safely omit such 
resonances. This justifies our neglect of the frame precession in the above calculation. 

However, the omission of the obliquity may have a serious effect on the results. Thus, we shall 
generalize Eqs. (38) so as to include the effect of the solar inclination and node in a Martian- 
centric frame. This generalized model based on the celebrated works by Kozai (1959) and Cook 
(1962) gives us: 



da 
Itt 
de 

'di 

di 

di 
doj 
Hi 







ISn'^eVl - e2 



An 



[2ABcos{2uj) - {A^ - B^) sin(2cj)] 



3n'^C 



4nVT 



—— = —Q cos i 



{ A [2 + 3e2 + Se^ cos(2cu)] + SEe^ sin(2cu) } 



(43a) 
(43b) 
(43c) 



Sn'^VT^ 
2n 



5AB sm{2u) + -{A^ - B^) cos(2cj) 



^ 3(^2 + B' 



+ 



15n a(A cos u + B sin cu) 



Anea' 



l-l{A^ + B') 



dQ 



3n"C 
An\J\ — sin; 



{SAe^ sin(2c^) + 5 [2 + 36^ - Se^ cos(2w)] } 



(43d) 
(43e) 
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where 



A 



= cos{n - Vl') cos(A') + cos(z') sin(A') sin(fi - Vl') 



(44a) 



B = cosi[- sm{n - n') cos(A') + cos(z') sin(A') cos{n - Q')] 



+ sin i sm(i') sm(A') 



(44b) 



C = sinz[cos A'sin(r2 — fi') — cos(z') sin(A') cos(r2 — f2') 



cos i sin i' sin A' 



(44c) 



Here, i' and Q' are the inchnation and right ascension of the ascending node of the Solar orbit 
in the Martian equatorial frame, respectively. The doubly averaged equations, obtained after 
averaging over the Sun's mean anomaly for a single period, are given by 



^ = 
dt 



(45a) 



dt An 



2v4Bcos(2cj) - (A2 _ 52) sin(2cj) 



(45b) 



di 



3n 



12 



dt Any/l^ 



{CA [2 + Se^ + 56^ cos(2cj)] + hCBe^ sin{2u)] 



(45c) 



duj 



-d cosi 



2n 



hABsm{2u) + -(^2 - 52) cos(2cu) - 1 



3(^2 + fi2j 



(45d) 



dVL 



3n 



12 



dt An\/l - e2 si 



sm« 



{5CAe2 sin(2cu) + CB[2 + Se^ - 56^ cos(2cu)] } 



(45e) 



the averaged quantities A^, B"^, AB, AC, BC being given by 
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= [s^,c^, - 0.5 s-,] + s^,sn'Sn'SnCn'Cf7 + 0.5 - 0.5s^,c^, (46a) 

= { [slc^, + 0.5s^,) - sn'S^,Sf7Cn'Cn + 0.58^,^, - 0.5 + 4} 

+ (cj'CnCn' + Ci/SnSnO Si'SiCi + 0.5s^, (46b) 

AB = {s^,sn'Cf7'C^ + [-s|sf7C^, + 0.5s^,sn] - 0.5s|sQ'Cn'} Ci 

+ 0.5ci/ ( sqCq' - cqSq') Sj/Si (46c) 

AC = 0.5 Ci> {sQCn' - cnsn') Sj'Cj 

+ {-s,^sn/ct7/c^ + [-s,^((sn)^c^,)] cq + 0.5s^,sn'Cn'} Sj (46d) 



2 



= (cj/cncn' + Cj/SnSQ') Si/Cj 
+ { [4 (cc)^ - 0-5 {si'f] - s^/Sf^/Sf^cn/Cf^ + -0.5s^,c^, - 4 + 0.5} s^Ci 
- 0.5 (cj/snsn' + Ci/CnCQ/) Si/ (46e) 

where we have used the compact notation C(.) = cos(-), S(.) = sin(-). To calculate the trigonometric 
functions of fl' and i' appearing in Eqs. (46), we shall utilise the geometry rendered Figure 1. 
The figure depicts the Martian spin axis, k, that is perpendicular to the equator of date, and the 
normal to the orbital plane (ecliptic of date), n. The Martian obliquity, e, is the angle between 
these two vectors, and is calculated based on the Colombo formalism: 

cos e = 'k-h = qSx+pSy + Fsz (47) 

where 



Sx = sin Ip sin hp , Sy = sin Ip cos hp , Sz = cos /p , F = \/l — — (48) 

and i' = e . 

The angle fi', lying in the equator of date, is subtended between the vectors LON12 and 
LON23. The first of these, LON12, points along the line-of- nodes obtained by the intersection 
of the equator of date and the invariable plane. This vector must be perpendicular to the plane 
defined by k and z (the normal to the invariable plane), so that 

LON21 = zxk=[sy, Sx , 0]"^ (49) 

The second vector, LON23 , is aimed along the line of nodes rendered by the intersection of the 
orbital plane and the equator of date. Based on the geometry of Fig. 1, we write: 

T 

LON23 = k X n = [ -SyF + SzP , SzQ - SxF , -SxP + Syq ] (50) 
By taking the scalar product of LON23 and LON21, we derive the direction cosine: 

^, LON23 • LON21 cos Jp cos i' - cos Jor-fe 

cosiZ = - — — ; — ; — — = 51 

ILON23I ILON21I sin Jp sin i' ^ ^ 

To derive an expression for sin Q' , we shall have to compute an auxiliary vector, j , which lies in 
the orbital plane and is normal to both k and LON21 (cf. Fig. 1): 

T 

j = k X LON21 = [ -SzSx , SzSy , 4 + 4 ] (^2) 
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The direction cosine between j and LON23 is then 



cos( ^-n') = LON.3-j^ ^53) 
^ 2 / ILON23I |j| 



Upon evaluating Eq. (53) we arrive at 

sinfi' = / ^"^'"^ (54) 



q cos hp — p sin hp 



sm^' 



Finally, the calculation of the orbital-elements' evolution is performed via Lagrange-type plane- 
tary equations, whose right-hand sides combine those of (1) - (5) and (45). 



2.5 The higher-order harmonics, the gravitational pull of the planets, 
the Yarkovsky effect, and the bodily tides 

In the current paper, we pursue a limited goal of taking into account the oblateness of the planet, 
its nonuniform equinoctial precession, and the gravitational pull of the Sun. These three items 
certainly do not exhaust the list of factors influencing the orbit evolution of a satellite. 

Among the factors that we intend to include into the model at the further stage of its de- 
velopment are the high-order zonal ( J3 , J4 , J| ) and tesseral ( C22 ) harmonics of the planet's 
gravity field, as well as the gravitational pull of the other planets, - factors whose role was com- 
prehensively discussed, for example, by Waz (2004). We also intend to include the bodily tides 
(Efroimsky & Lainey 2007) and the Yarkovsky effect (Nezvorny & Vokrouhlicky 2007) - factors 
whose importance increases at long time scales . 



3 Comparison of a purely numerical and 
a semianalytical treatment of the problem 

One of our goals is to check the applicability limits (both in terms of the initial conditions and 
the permissible time scales) of our semianalytical model written for the osculating elements in- 
troduced in a frame co-precessing with the equator of date. This check will be performed by an 
independent, purely numerical, computation that will be free from whatever simplifying assump- 
tions (all terms kept, no averaging performed.) The straightforward simulation will be carried out 
in terms of Cartesian coordinates and velocities defined in an inertial frame of reference. Both 
the semianalytical calculation of the elements in a co-precessing frame and the straightforward 
numerical integration in inertial Cartesian axes will be carried out for Deimos. 



3.1 Integration by a purely numerical approach. 

The numerical integration of Deimos' orbit can be performed using Cartesian coordinates defined 
relatively to the Solar system invariable plane. As we also have to compute the Martian polar 
axis motion, there are two vector differential equations to integrate simultaneously. One is the 
Newton gravity law: 
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where Vf/ has components 



dM 



dyU 



pi -h 




ri5 








PIJ2 


-( 


as 




T \ 




pI J2 


Z 






r 





sm 



sm 



sm 



3 
2 

3 
2 

3 
2 



3 sin I 



sin/p sin hp 



3 sin0 sin/p cos hp 



3 sin 6 cos In 



(56) 



Here and r = [x, y, z) denote, correspondingly, the latitude of Deimos relative to the Martian 
equator and the position vector of Deimos related to the Martian center of mass; pe is the Martian 
equatorial radius; M and m stand for the masses of Mars and Deimos, respectively. Angles hp 
and Ip are the longitude of the node and the inclination of the planet's equator of date relative 
to the invariable plane (see Section 2.2). Integration in this, inertial, frame offers the obvious 
advantage of nullifying the inertial forces. 

Table 2 gives the initial conditions for our simulation, expressed in terms of the Keplerian 
orbital elements. Table 3 presents these initial conditions in a more practical form, i.e., in terms 
of the Cartesian positions and velocities corresponding to the said elements. A transition from 
the Keplerian elements to these Cartesian positions and velocities is a two-step process. First, 
we take orbital elements defined in a frame associated with the Martian equator of date (i.e., in 
a frame co-precessing but not co-rotating with the planet) and transform them into Cartesian 
coordinates and velocities defined in that same frame. Then, by two successive rotations of angles 
— Ip and — hp , we transform them into Cartesian coordinates and velocities related to the 
invariable plane. These initial positions and velocities were used to begin the integration. 

At each step of integration of (55), the same two rotations are performed on the components 
Vf/ given by (56). As mentioned above, to afford the absence of inertial forces on the right-hand 
side of (55) one must write down and integrate (55) in the inertial frame. Since the analytical 
expressions (56) for VU contain the latitude , they are valid in the co-precessing coordinate 
system and, therefore, need to be transformed to the inertial frame at each step. To carry out 
the transformation, one needs to know, at each step, the relative orientation of the Martian polar 
axis and the inertial coordinate system. The orientation is given by the afore mentioned Colombo 
model. This is how our second equation, the one of Colombo, comes into play: 



^ = a(k ■ n)(k X n) 

dt V A ; 



(57) 



All in all, we have to integrate the system (55 - 57). Table 4 gives the initial conditions used for 
integrating (57), while Table 5 gives the numerical values used for the parameters involved. It is 
worth noting that the values of Table 4 were calculated based on Ward (1974). 

The software used for numerical integration of the system (55 - 57) is called NOE (Numerical 
Orbit and Ephemerides) , and is largely based on the ideas and methods developed in Lainey, 
Duriez & Vienne (2004). This numerical tool was created at the Royal Observatory of Belgium 
mainly for computations of the natural satellite ephemerides. It is an A^-body code, which in- 
corporates highly sensitive modelling and can generate partial derivatives. The latter are needed 
when one wants to fit the initial positions, velocities, and other parameters to the observation 
data. To save the computer time, an optimised force subroutine was built into the code, specif- 
ically for integrating the above equations. This appliance, based on the RA15 integrator offered 
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by Everhart (1985), was chosen for its speed and accuracy. During the integration, a variable 
step size with an initial value of 0.04 day was used. To control the numerical error, back and 
forth integrations were performed. In particular, we carried out a trial simulation consisting of a 
thousand-year forward and a subsequent thousand-year back integration. The satellite displace- 
ment due to the error accumulated through this trial was constrained to 150 meters. Most of this 
150 meter difference comes from a numerical drift of the longitude, while the numerical errors in 
the computation of the semi-major axis, the eccentricity, and the inclination were much lower. 
These errors were reduced for this trial simulation to only 10~^ km , 10"^" , and 10"^'' degree, 
respectively. This provided us with a high confidence in our subsequent numerical results. 

As a complement to the said back-and-forth check, the energy-conservation criterium was 
used to deduce, in the first approximation, an optimal initial step-size value and to figure out the 
numerical error proliferation. (It is for this energy-conservation test that we introduced a non- 
zero mass for Deimos. Its value was taken from Smith, Lemoine & Zuber (1995).) Applicability 
of this criterium is justified by the fact that the numerical errors are induced mostly by the fast 
orbital motion of the satellite.^ 



parameters 


numerical values 


a 


23459 km 


e 


0.0005 


i 


0.5 deg and 89 deg 




10 deg 




5 deg 


M 


deg. 



Table 2: The orbital elements values taken as initial conditions for our simulations. 



satellite 




X 


y 


z 


position km {i 


= 0.5) 


22648.3376439 


6068.52353055 


17.8332361962 


velocity km/s i 


[i = 0.5) 


-0.349882011871 


1.30576017694 


1.175229063323 x 10"^ 


position km {i 


= 89) 


22996.9921622 


4091.20549954 


2043.25303109 


velocity km/s i 


[i = 89) 


-0.120115009144 


2.686751629968 x lO^^ 


1.34652528539 



Table 3: The initial positions and velocities used for Deimos. The first two rows correspond to 
the low-inclination case (0.5 degree); the last two rows correspond to the high-inclination case 
(89 degrees). 



^Although the planet-satellite system is subject to an external influence (the solar torque acting on the planet), 
over short time scales this system can be assumed closed. In order to check the integrator efficiency and to determine 
an optimal initial step size, we carried out auxiliary integrations of (55) - (56), with the Colombo equation (57) 
neglected and with the energy presumed to conserve. These sevcral-thousand-year-long trial integrations, with the 
energy-conservation criterium applied, led us to the conclusion that our integrator remained steady over long time 
scales and that the initial step of 0.04 day was optimal. Then this initial step size was employed in our integration 
of the full system (55) - (57). 
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parameters 


numerical values 




25.25797549 deg 


hpito) 


332.6841708 deg 



Table 4: Initial conditions used for the integration of Eq. (57). These values were calculated 
based on Ward (1974). 



parameters 


numerical values 


Martian mass (GM) 


42830 km^s"^ 




1960.45 X 10-6 


Equatorial radius 


3397 km 


Deimos mass 


0.091 X 10^3 km^s-^ 


a 


3.9735 X 10"^ rad/yr 



Table 5: Parameter values used in our simulations. 



3.2 Integration by the semianalytical approach 

The theory of satellite-orbit evolution, based on the planetary equations whose right-hand sides 
combine those of (1 - 5) and (45), is semianalytical. This means that these equations for the 
elements' secular parts are derived analytically, but their integration is to be performed numeri- 
cally. This integration was carried out using an 8*'*-order Runge-Kutta scheme with relative and 
absolute tolerances of 10"^^ . Kilograms, years, and kilometers were taken as the mass, time, 
and length units, correspondingly. 



3.2.1 Technicalities 



To integrate the planetary equations, one should know, at each time step, the values of Ip and 
hp , which are the time derivatives of the inclination and of the longitude of the node of the 
equator of date with respect to the equator of epoch. These derivatives will be rendered by the 
Colombo equation (17), after formulae (15 - 16) get inserted therein: 



a 



( sin Ip sin hp cos hp — q p sin Ip + 2 p q sin Ip cos^ hp 



— p^ sin Ip cos hp sin hp + q \/\ — q^ — p^ cos Ip cos hp 



V 



and 
hr, = 



a 



\/\ — (f' — p^ cos Ip sin hp j , 

( — g -|- 2g cos^ /p) cos^ hp — 2q cos^ Ip + q 



(58) 



{p — 2p cos^ Jp) cos hp 



sin hr 



p2 _ q2 



sin Ir, 



2 qp 



p^) cos Jp cos^ /ip + { — p' 

cos Ip cos^ hp — 2pq cos 
sin hp 



2q^ + l) cos Jp 
ihp cos Ip I 



(59) 
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Equations (58) and (59) are then integrated (with the initial conditions /p(to) and hpito) bor- 
rowed from Table 4) simultaneously with the planetary equations (1-5). Through formulae (6), 
the above expressions for Ip and hp yield the expressions for the components of . As 
can be seen from (12 - 14), integration of the planetary equations also requires the knowledge 
of the derivatives /ii , /i2 , and /is at each integration step. These can be readily obtained by 
differentiating (6). The resulting closed- formed expressions for fii , fi2 , jj^z are listed in Appendix 
A. The final step required for numerical integration of the planetary equations is substitution of 
formulae (9 - 12), with the initial conditions from Table 2. 

3.2.2 The plots and their interpretation 

Fig. 2 depicts the history of the planetary equator, in the Colombo approximation, over 1 Byr. 
The inclination exhibits long-periodic oscillations bounded within the range of 20.3 deg < Ip < 
30.3 deg , while the node regresses at a rate of hp = 0.00202 deg /yr . The obliquity, e, varies in 
the range 15.2 deg < e < 35.5 deg . A magnified view of the obliquity for 5 Myr is shown in 
Fig. 3. The time history of the obliquity closely matches the results reported by Ward (1974). 

Integration of our semianalytical model gives plots depicted in Figures 4 and 5, for a low 
initial inclination, and in Figures 6-7, for a high initial inclination. From Fig. 4 we see that the 
variable equinoctial precession does not inflict considerable changes upon the satellite's inclination 
relative to the precessing equator of date. The orbit inclination remains bounded within the region 
0.3 deg < i < 2.5 deg . This means that the "Goldreich lock" (inclination "locking," predicted 
by the Goldreich (1965) model for small inclinations and for uniform equinoctial precession) works 
also for nonuniform Colombo precession of the equator. 

However, this observation no longer holds for large initial inclinations. In the iq = 89 deg case, 
it is clearly seen that the node is greatly affected by the presence of the equinoctial precession. The 
equinoctial precession also affects the magnitude of the inclination variations. This case exhibits 
chaotic dynamics that are sensitive to any additional perturbing inputs. Fig. 6 shows that without 
precession the inclination of Deimos' orbit varies within the range of 84.5 deg < z < 95 deg . 
With the precession included, the inclination gains about one degree in amplitude, varying in the 
range 83.5 deg < i < 96 deg . This effect is accentuated when examining a magnification of the 
inclination over a 5 Myr span, as shown in Fig. 6. The chaotic nature of the inclination is clearly 
seen. The irregular dynamics is characterised by chaotic switches between the maximum and 
minimum inclination values, a phenomenon referred to as "crankshaft" , to be further discussed 
in the sequel. 

Both in the near-equatorial case (as in Fig. 5) and the near-polar case (as in Fig. 7), variations 
of the semimajor axis are of order 10~^ % . This smallness is in compliance with formula (44) 
in Efroimsky (2006a,b), according to which the changes of a generated by the variations of the 
equinoctial precession rate are extremely small. (Formula (38a) from subsection 2.4 above tells 
us that the direct pull of the Sun exerts no influence upon a at all.) 

Similarly, in both cases (Fig. 5 and Fig. 7) the variations of eccentricity, remain small, about 
10"^ % . While formula (44) in Efroimsky (2006a,b) promises to e only tiny variations due to 
precession, our formula (38b) from subsection 2.4 above gives to e a slightly higher variation 
rate, rate that still remains insufficient to raise the quasiperiodic changes in e above a fraction 
of percent. 

The plots in Figs. 5 and 7 depict also the time evolution of uj . The line of apsides steadily 
regresses in the near-polar case and steadily advances in the near-equatorial case. 

To examine the precision of our semianalytical model, we have compared the results of its 
integration with the results stemming from a purely numerical simulation performed in terms of 
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inertial Cartesian coordinates and velocities (see subsection 3.1). The necessity for this check 
was dictated, mainly, by the fact that within the semianalytical model the short-period terms are 
averaged out,^ while the straightforward numerical integration of (55 - 56) neglects nothing. We 
carried out comparison of the two methods over 10 Myr only. The outcomes, both for = 0.5 
degrees and = 89 deg , were in a good agreement. As an example, the top plot in Fig. 8 
shows the comparison of the inclination evolution calculated by the semianalytical and purely 
numerical methods over 10 Myr in the case of = 0.5 deg . The bottom plot in Fig. 8 depicts a 
similar comparison for the case of ^o = 89 deg. Since the inclination exhibits chaotic behaviour, 
comparing the semianalytical and purely numerical calculations point-by-point (i. e., computing 
the differences between these two signals) would not be useful. Therefore, we chose to compare 
the mean and standard deviation (STD) of the inclination, in the semianalytical and purely 
numerical simulations. We also compared the extremum values. The results of this comparison 
are summarised by Table 6. The table quantifies the agreement between the models depicted by 
Fig. 8. It clearly demonstrates that the two simulations agree up to fractions of a percent. 

All in all, the outcome of our computations is two-fold. First, we have made sure that the 
semianalytical model perfectly describes the dynamics over time scales of, at least, dozens of 
millions of years. Stated differently, the short-period terms and the terms of order 0(/x^) play no 
role over these time spans. Second, we have made sure that the "Goldreich lock" initially derived 
for very low inclinations and for uniform equinoctial precession, works well also for variable 
precession, though for low initial inclinations only. 



«o [deg] 


Model 


STD [deg] 


Mean [deg] 


Max. Value [deg] 


Min. Value [deg] 


0.5 


Semianalytical 


0.60 


1.519 


2.45 


0.3063 




Cartesian 


0.601 


1.53 


2.465 


0.3056 


89 


Semianalytical 


3.10 


90.085 


95.9713 


84.027 




Cartesian 


3.09 


89.92 


95.9769 


84.0054 



Table 6: Statistical properties of the inclination. Comparison of the results of the semianalytical 
and the purely numerical computation. 



3.2.3 Looking for trouble 

A natural question arises as to whether the considered examples are representative. One may 
enquire if, perhaps, there still exists a combination of the initial conditions yielding noticeable 
variations of the satellite orbit inclination during the primary's variable equinoctial precession over 
vast spans of time. To answer this question, we should scan through all the possible combinations 
of initial conditions, to identify a particular combination that would entail a maximal inclination 
excursion relative to the initial inclination. Stated more formally, we should seek a set of initial 
conditions { h*Q , /*q , iQ , Qq, Uq} maximising the objective function \i{t) — zq | : 

{^pO' ^pO' ^0' ^0. ^o) = argmax|i(t) - io\ (60) 

This problem belongs to the realm of optimisation theory. The optimisation space is constituted 
by the entire multitude of the permissible initial conditions. The sought after combination of 

^ We remind that in equations (1 - 5) the exact /x-dependent terms are substituted with their orbital averages 
(12 - 14). 
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initial conditions will be called the optimal set. 

In this situation, the traditional optimisation schemes (such as the gradient search or the 
simplex method) may fail due to the rich dynamical structure of our problem - these methods 
may lead us to a local extremum only. Thus, the search needs to be global. It can be carried 
out by means of Genetic Algorithms (GA's) [20, 22, 23, 24]. These are wont to supersede the 
traditional optimisation procedures in the following aspects. First, instead of directly dealing 
with the parameters, the GA's employ codings (usually, binary) of the parameter set ("strings," 
in the GA terminology). Second, instead of addressing a single point of the optimisation space, 
the GAs perform a search inside a population of the initial conditions. Third, instead of pro- 
cessing derivatives or whatever other auxiliary information, the GAs use only objective-function 
evaluations ("fitness evaluations"). Fourth, instead of deterministic rules to reiterate, the GAs 
rely upon probabilistic transition rules. Additional details on the particular GA mechanism used 
herein can be found in Appendix B. 

A GA optimisation was implemented using the parameter values given in Table 7. The search 



parameters 


numerical values 


Population size 


30 


Number of generations 


100 


String length 


16 bit 


Probability of crossover 


0.99 


Probability of mutation 


0.02 



Table 7: Parameter values used for the GA optimization 

for the inclination-maximising initial conditions resulted in the following set: 

/;o = 72.5deg, /i^o = 211.324 deg, ^o = 100.543 deg , = 111-538 deg , cj* = 234.913 deg . 

Thus, the initial orbit is retrograde and, not surprisingly, near-polar. The resulting time histories 
for a 0.2 Byr integration are depicted in Fig. 9, for i and Q . In both cases the inclination 
amplitude is relatively large: The inclination varies within the range of 79 deg < z < 102 deg. 

This example clearly shows that the equinoctial precession is an important effect for evolution 
of satellite orbits. As shown in the upper pane of Fig. 9, had we neglected the precession, the 
magnitude of the oscillations would be about twice smaller. The inclusion of the precession in 
the model qualitatively modifies the behavior, inducing large-magnitude chaotic variations of the 
inclination, a phenomenon that cannot be detected without including the precession alongside 
the oblateness and solar gravity. 

4 Comparison of the semianalytical results 
with those rendered by Goldreich's model 

The final step in our study will be to compare the semianalytical model to Goldreich's approx- 
imation (22 - 26). To that end, we integrate our semianalytical model for 20 Myr, using the 
initial conditions from Table 2 with = 89 deg; and compare the outcome with that resulting 
from Goldreich's approximation simulated with the same initial conditions. The results of this 
comparison are depicted in Figs. 10 - 11. Specifically, Fig. 10 compares the time histories of 
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Ip and hp . There are noticeable differences in the dynamics of Ip . While Goldreich's approx- 
imation assumes a constant Ip , the semianalytical model is based on the Colombo calculation 
of the equinoctial precession, calculation that predicts considerable oscillations within the range 
21 deg < /p < 30 deg . Beside this, in our semianalytical model we take into account the direct 
gravitational pull exerted by the Sun on the satellite. All this entails differences between the 
dynamics predicted by our semianalytical model and the dynamics stemming from the Goldreich 
approximation. These differences, for i and fl , are depicted in Fig. 11. In Goldreich's model, 
i stays very close to the initial value: 88.27 deg < i < 89.01 deg , a behaviour that makes 
the essence of the Goldreich lock. However, in the more accurate, semianalytical model we have 
84 deg < « < 96 deg . The time history of Q , too, reveals that Goldreich's approximation does 
not adequately model the actual dynamics, since it predicts a much larger secular change than 
the semi-analytical model. 

All in all, the dynamics (i.e., particular trajectories) predicted by the two models are quanti- 
tatively different. At the same time, when it comes to the most physically important conclusion 
from the Goldreich approximation, the "Goldreich lock" of the inclination, one may still say that, 
qualitatively, the semianalytical model confirms the locking even in the case when the equinoctial 
precession is variable, and the solar pull on the satellite is included. The locking survives even for 
highly inclined satellite orbits. It is, though, not as stiff as predicted by the Goldreich model: we 
can see from Fig. 11 that the orbit inclination varies within a five-degree span, while the Goldreich 
approximation would constrain it to fractions of a degree. 

An intriguing fine feature of the inclination evolution, which manifests itself for orbits close to 
polar, is the "crankshaft". This kind of behaviour, well defined in Fig. 11, is not rendered by the 
Goldreich model, because that model was initially developed for small inclinations and uniform 
precession. One might suspect that the "crankshaft" is merely a numerical artefact, had it not 
been discovered under different circumstances (in the absence of precession but in the presence 
of a third body) by Zhang & Hamilton (2005). This kind of pattern may be generic for the close 
vicinity of z = 90 deg . 

5 Conclusions 

In the article thus far, we continued developing a tool for exploring long-term evolution of a 
satellite orbit about a precessing oblate primary. In particular, we were interested in the time- 
dependence of the orbit inclination relative to the moving equator of date. Our model includes 
three factors: J2 of the planet, the planet's nonuniform equinoctial precession described by 
the Colombo formalism, and the gravitational pull exerted by the Sun on the satellite. The 
problem was approached using different methods. One, semianalytical, was based on numerical 
integration of the averaged Lagrange-type equations for the secular parts of the Keplerian orbital 
elements introduced in a noninertial reference frame coprecessing with the planetary equator of 
date. The right-hand sides of these equations consisted of precession-generated contributions and 
contributions due to the direct pull of the Sun. The other approach was a straightforward, purely 
numerical, computation of the satellite dynamics using Cartesian coordinates in a quasi-inertial 
reference frame. 

The results of both computations have demonstrated a good agreement over 10 Myr. This 
means that the semianalytical model adequately describes the dynamics over this time span. 
Specifically, the terms neglected in the semianalytical model (the short-period terms and the 
terms of order 0(/x^) ) play no significant role on this time scale. 

Our calculations have shown the advantages and the limitations of a simple model developed 
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by Goldreich (1965) for uniform equinoctial precession and low inclinations. Though his model 
does not adequately describe the dynamics (that turns out to be chaotic), the main physical 
prediction of Goldreich's model - the "Goldreich lock" - sustains the presence of the Sun and 
variations of equinoctial precession, provided the initial inclination is sufficiently low. For low 
initial inclinations, the inclination exhibits variations of order fractions of a degree. For higher 
inclinations, however, it varies already within a span of about ten degrees. For near-polar orbits, 
the inclination behaviour demonstrates the "crankshaft", an chaotic pattern not accounted for 
by the Goldreich model. The "crankshaft" emerges because in our model both the precession 
variations and the pull of the Sun are included into the model. However, numerical experiments 
have also shown that the "crankshaft" gets generated by each of these two factors separately. 



Appendix A: Closed-Form Expressions for /ii, /i2, /is 

Using the compact notation C(.) = cos(-) and S(.) = sin(-) , and conforming to the procedure 
described in the text, we obtain the following expressions for fii, fi2, As : 
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and 
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Appendix B: Niching Genetic Algorithms 

The most commonly used Genetic Algorithm (GA) is the so-called "Simple GA." To perform an 
evolutionary search, the Simple GA uses the operators of crossover, reproduction, and mutation. 
A crossover is used to create new solution strings ("children" or "offspring") from the existing 
strings ("parents"). Reproduction copies individual strings according to the objective function 
values. Mutation is an occasional random alteration of the value of a string position, used to 
promote diversity of solutions. 

Although Simple GA's are capable of detecting the global optimum, they suffer from two 
main drawbacks. First, convergence to a local optimum is possible due to the effect of premature 
convergence, where all individuals in a population become nearly identical before the optima has 
been located. Second, convergence to a single optimum does not reveal other optima, which may 
exhibit attractive features. To overcome these problems, modifications of Simple GA's were con- 
sidered. These modifications are called niching methods, and are aimed at promoting a diversity 
of solutions for multi-modal optimisation problems. In other words, instead of converging to 
a single (possibly local) optimum, niching allows for a number of optimal solutions to co-exist, 
and it lets the designer choose the appropriate one. The niching method used throughout this 
study is that of Deterministic Crowding. According to this method, individuals are first randomly 
grouped into parent pairs. Each pair generates two children by application of the standard genetic 
operators. Every child then competes against one of his parents. The winner of the competition 
moves on to the next generation. By using the notation Pj for a parent, Cj for a child, /(■) for 
a fitness, and d{-) for a distance, a pseudo-code for the two possible parent-child tournaments 
can be written as follows: 



If [d(Pi, Ci) + t/(P2, C2) = d{Pi, C2) + d{P2, Ci)] 

If /(Ci) > /(Pi) replace Pi with Ci 
If /(C2) > /(P2) replace P2 with C2 
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Else 



If /(Ci) > /(P2) replace P2 with Ci 
If /(C2) > /(Pi) replace Pi with C2 

In addition to applying the Deterministic Crowding niching method, we used a two-point 
crossover instead of a single-point one. In the Simple GA, the crossover operator breaks the 
binary string of parameters, the "chromosome," at a random point and exchanges the two pieces 
to create a new "chromosome." In a two-point crossover, the "chromosome" is represented with 
a ring. The string between the two-crossover points is then exchanged. The two-point crossover 
or other multiple-point crossover schemes have preferable properties when optimisation highly 
nonlinear functions is performed. 
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Figures 




Figure 1: The geometry of the precessing Mars. The Martian spin axis, k, is perpendicular to the 
equator of date, and the normal to the orbital plane (ecliptic of date), n. The Martian obliquity, 
e, is the angle between these two vectors. 
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Figure 2: Evolution of the Martian inclination, the longitude of the node of the equator of date 
relative to that of epoch, and the obliquity over 1 Byr, obtained in the Colombo approximation. 
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Figure 3: Evolution of the the Martian obhquity for 5 Myr, calculated through formula (47). This 
curve closely matches the result of Ward (1974). 
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Figure 4: Evolution of the inclination (initially set to 0.5 degree) and of the longitude of the node 
of Deimos over 1 Byr. The plot, obtained by integration of the semianalytical model, exhibits 
inclination locking and a uniform regression of the node. 
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Figure 5: Evolution of the semimajor axis, eccentricity, and argument of periapsis of Deimos over 
1 Byr. (The inchnation was initially set to 0.5 degree.) Both the semimajor axis and eccentricity 
exhibit quasiperiodic motion about their initial values. (The variations of the semimajor axis and 
eccentricity are so small that it is more convenient to plot a — oq and e — Cq .) The plots were 
obtained by integration of the semianalytical model. 
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Figure 6: Evolution of the inclination and of the longitude of the node of Deimos over 1 Byr. (The 
inclination was initially set to 89 degrees.) The plot, obtained by integration of the semianalytical 
model, exhibits inclination variation in the range ±5 deg and a chaotic evolution of the node; if 
precession is neglected, the inclination oscillations are smaller in magnitude. The chaotic nature 
of the inclination variation is referred to as "crankshaft chaos" . 
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Figure 7: Evolution of the semi-major axis, eccentricity, and argument of periapsis of Deimos 
over 1 Byr. (The inchnation was initially set to 89 degrees.) The semi-major axis exhibits 
chaotic behavior and the eccentricity exhibits long-periodic motion about the initial value. (The 
variations of the semimajor axis are so small that it is more convenient to plot a — Qo rather 
than a .) The plots were obtained by integration of the semianalytical model. 



35 




Cartesian 

Semianalytical 




2 4 6 8 10 

Time[yr] 



Figure 8: Comparison of the semi- analytical model to a purely numerical integration. The plot 
shows the inclination as predicted by the two models. The top plot shows the case of the initial 
inclination zq = 0.5 deg. In this case, the standard deviation of the inclination in the semian- 
alytical and inertial models differs by 0.175%, and the mean value by 0.77%. The bottom plot 
shows the case of the initial inclination io = 89 deg. In this case, the differences between the 
semianalytical and the numerical model amount to 0.4% in standard deviation and 0.18% in the 
mean value. 
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Figure 9: Evolution of the inclination and longitude of the node, for the initial conditions that 
entail maximal variations of the orbit inclination. The plot, obtained by integration of the semi- 
analytical model, demonstrates that the inclination, initially set at 100.5 deg, plunges to less than 
80 deg and then returns to its initial value. The switches between maximum and minimum values 
exhibit "crankshaft" chaos. Without precession, the inclination magnitude is much smaller, and 
the node regresses uniformly. 
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Figure 10: Preparation to comparison of the semianalytical model to Goldreich's model. In 
Goldreich's model, the Martian equator is assumed to precess uniformly, while the semianalytical 
model takes into account (via Colombo's equation) variations of the equinoctial precession. 
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Figure 11: Comparison of the semi-analytical model to Goldreich's model for a satellite of Deimos' 
mass and with the initial conditions given by (61). While in Goldreich's model the inclination 
remains tightly locked, the semianalytical model reveals much larger variations of i . 
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